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ABSTRACT 


Two-phase  Hartmann  flows  in  an  MHD  generator  duct  of  rectangular  cross- 
section  are  examined  and  numerical  values  for  the  velocity  fields  and  induced 
current  streamlines  for  fluids  of  spatially-dependent  electrical  conductivity 
are  determined  as  a function  of  load  factor,  Hartmann  number  and  channel  aspect 
ratio.  This  study  is  of  considerable  practical  interest  in  connection  with  some 
experimental  work  in  progress  at  Argonne  Laboratories  on  two  phase  liquid  metal- 
gas  MHD  generators  and  the  results  presented  herein  may  enhance  the  understanding 
of  the  operating  characteristics  of  such  devices.  The  problems  considered 
include:  1)  one-dimensional  Hartmann  flow  between  two  parallel  plates  at  high 

Hartmann  number  when  the  fluid  electrical  conductivity  decreases  with  increasing 
distance  from  the  walls,  2)  two-dimensional  Hartmann  flow  in  rectangular 
channels  of  finite  aspect  ratio  having  two  insulating  walls  and  two  thin  walls 
of  finite  conductivity  representing  the  electrodes  and,  3)  compressible  one- 
dimensional isothermal  MHD  flow  with  zero-slip  between  its  liquid  and  gas 
components.  Results  indicate  that  conductivity  gradients  can  produce  inflected 
velocity  profiles  which  may  be  hydrodynamically  unstable,  that  for  larger  load 
factors  the  shunt  current  in  the  Hartmann  boundary  layers  become  considerable 
and  that  two-phase  MHD  flows  of  moderate  void  fractions  have  sound  speeds 
considerably  lower  than  in  its  components.  Turbulent  effects  are  not  considered 
in  these  studies  since  turbulent  fluctuations  are  expected  to  be  suppressed 
for  the  high  Hartmann  number  flows  of  interest  in  two-phase  MHD  generators. 
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I . INTRODUCTION 

An  analytical  study  on  two-phase  Hartmann  flows  in  the  MHD  generator 
configuration  has  been  in  progress  here  at  the  University  of  Florida  since 
December  1975.  The  objective  of  this  ONR  sponsored  program  has  been  to 
better  understand  the  behaviour  of  electrically  conducting  flows  with  spatially 
dependent  electrical  conductivity  in  the  presence  of  externally  applied 
magnetic  fields  and  in  particular  to  aid  in  the  understanding  of  some  of  the 
experimental  observations  which  have  been  and  are  being  obtained  in  a major 
experimental  research  program  at  Argonne  Laboratories  on  liquid  metal-gas 
MHD  power  generation [ 1 , 2).  The  hydromagnetic  flow  behaviour  in  such  two- 
phase  devices  is  of  special  interest  to  the  Navy  as  these  generators,  if 
they  can  be  made  to  function  in  scaled-up  versions,  have  good  efficiencies 
and  high  power  densities  which  would  make  them  ideally  suited  for  on-board 
applications.  At  the  present  time  an  effort  is  underway  at  Argonne 
Laboratories  to  construct  a two-phase,  sodium-argon  MHD  generator  with  a net 
30  MW  power  output [1}.  It  is  part  of  our  effort  at  the  University  of 
Florida  to  gain  a better  understanding  of  the  properties  of  the  MHD  flows 
encountered  in  such  generators  and  possibly  to  aid  in  the  surmounting  of 
expected  difficulties  introduced  by  upscaling,  end  effects  and  liquid-gas 
separation  within  the  generator  channel.  We  have  made  an  effort  to  initiate 
and  maintain  communications  between  the  various  research  groups  working  in 
liquid  metal-gas  MHD  power  generation.  Discussions  to  date  have  been  held, 
among  others, with  M.  Petrick's  Group  at  Argonne  Laboratories  and  with  H. 
Branover  of  Ben-Gurion  University. 

The  present  report  gives  a summary  of  results  obtained  by  us  here 
at  the  University  of  Florida  during  the  first  year  of  the  contract  period. 
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The  personnel  Involved  In  the  combined  analytical-numerical  research 
effort  have  been  Drs.  E.  R.  Lindgren  and  U.  H.  Kurzweg  as  co-principal 
investigators  , together  with  R.  E.  Elkins  as  a full-time  post-doctoral 
investigator  and  T.  A.  Trovillion  as  a PhD  canidate.  The  initial  phase 
of  our  study  included  a careful  examination  of  the  relevant  literature  in 
the  field  and  a formulation  of  an  appropriate  set  of  equations  governing 
two-phase  magnetohydrodynamic  flows.  This  was  followed  by  major  analytical 
efforts  in  the  areas  of:  1)  One-dimensional  Hartmann  flow  for  fluids  of 

variable  electrical  conductivity,  2)  Two-dimensional  Hartmann  flow  in 
channels  of  rectangular  cross-section  and  3)  One-dimensional  compressible 
MHD  flow  under  zero-slip  conditions.  Details  of  our  results  in  these 
areas  are  presented  below. 
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II.  ONE-DIMENSIONAL  HARTMANN  FLOW  WITH  VARIABLE  ELECTRICAL  CONDUCTIVITY 


L 


Our  initial  studies  involved  the  determination  of  velocity  and  current 
distributions  in  one-dimensional  Hartmann  flow  between  two  parallel  plates 
when  the  fluid  electrical  conductivity  and  viscosity  was  a function  of  the 
normal  distance  from  the  bounding  walls  and  a constant  external  magnetic 
field  was  applied.  This  represents  a first  approximation  to  homogeneous, 
laminar,  two-phase  flow  in  an  MHD  channel  of  large  aspect  ratio  when  bubble 
motion  relative  to  that  of  the  liquid  is  neglectad.  The  appropriate 
differential  equation  for  such  a laminar,  one-dimensional  flow  in  this 
geometry  can  be  derived  from  the  standard  continuum  MHD  equations  as  given 
for  example  by  Hughes  and  Young  [3]  and  reads 
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where  p(y)  is  the  viscosity  and  a(y)  the  electrical  conductivity  both  of 
which  are  taken  as  functions  of  the  normal  distance  from  the  walls.  Here  z 
is  the  flow  direction  coordinate,  W(y)  the  axial  velocity,  the  constant 
applied  magnetic  intensity  in  the  y direction,  E the  electric  field  in 
the  current  flow  direction  x and  3p/3z  the  constant  pressure  gradient.  As 
we  were  attempting  to  approximate  the  experimental  conditions  encountered 
in  the  Argonne  Group's  studies  [2],  the  electrical  conductivity  was  assumed 


to  have  the  void  related  y dependence  given  by  the  empirical  formula 
a = oq  [0.97083  - 1.04679  a - 2.22367a2  + 3.52146a3  - 1.11742a4] 


(2) 


with  a parabolic  void  distribution 

a = 0.7  - 0.5(y/H)2  = Gas  Volume/Total  Volume  (3) 

corresponding  to  that  found  at  station  4 in  their  non-tapered  channel 
experiments.  Here  y=+H  are  the  locations  of  the  bounding  walls.  Note  again 
that  in  this  model  the  two-phase  liquid  metal-gas  mixture  is  being  treated 
as  a homogeneous  fluid  with  variable 
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conductivity  properties.  Although  such  an  approach  is  certainly  not  exact, 
it  nevertheless  can  be  expected  to  shed  valuable  new  light  on  the  problem. 

Also  it  should  be  noted  that  we  are  applying  Eq.  (2)  to  yield  a spatial  vari- 
ation for  o * while  the  original  Argonne  Laboratory  results  [2]  gave  this 
conductivity  behaviour  as  an  averaged  value  across  the  channel  in  the 
absence  of  a magnetic  field.  Because  of  this  difference,  we  may  treat 
Eq.(2)  as  one  among  many  experimental  conduction  formulas  showing  a 
monotonically  decreasing  electrical  conductivity  with  increasing  distance 
from  the  bounding  walls  of  the  channel.  The  viscosity  coefficient  in  our 
calculations  is  taken  as  constant.  This  simplifying  assumption  is  reasonable 
as  viscous  effects  can  be  expected  to  be  confined  to  a very  narrow  region 
near  the  walls  for  the  high  Hartmann  number  cases  of  interest.  In  such  a 
thin  Hartmann  layer  the  viscosity  coefficient  can  be  treated  essentially  as 
constant  although  it  may  vary  considerably  in  the  core  of  the  flow.  The 
boundary  conditions  for  the  problem  are  the  no  slip  velocity  condition  W(+H)=0 
and  the  integral  requirement 
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where  R is  the  external  applied  resistance,  L is  the  channel  length  in  the 
flow  direction  and  2b  is  the  distance  between  electrodes  in  the  x direction. 

Equation  (1),  subjected  to  the  stated  boundary  conditions,  was 
solved  by  both  an  analogue  and  digital  numerical  approach  since  analytic 
solutions  were  not  found  to  be  possible  for  the  conductivity  distribution 
given  by  Eqs.  (2)  and  (3)  except  within  the  core  of  the  flow.  The  analogue 
solutions  were  obtained  for  Hartmann  numbers  M=BCH  J o 0 / u " up  to  10  but  the 
solution  became  unstable  for  still  larger  values.  This  difficulty  is  due 
to  the  rapid  change  in  magnitude  of  the  velocity  profile  in  the  wall  near 
region  for  larger  M and  led  us  to  abandon  the  analogue  calculations  in 
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favor  of  the  finite  difference  digital  approach.  The  digital  calculations 


were  carried  out  by  dividing  the  channel  half-width  0<y<H  into  ten  equally 
spaced  subintervals  within  the  Hartmann  boundary  layer  located  approximately 
in  H>y>H[ 1- (1/M) ] and  ten  equally  spaced  subintervals  between  the  boundary 
layer  edge  and  the  channel  axis  at  y=0.  In  this  manner  a simple  matrix 
inversion  yielded  solutions  for  W(y)  for  any  desired  Hartmann  number  and 
calculations  for  both  constant  conductivity  and  the  variable  distribution 
given  by  Eqs.  (2)  and  (3)  were  made  for  Hartmann  numbers  as  high  as  M=1000 
without  encountering  convergence  difficulties.  Typical  results  for  the 
velocity  variation  in  the  variable  conductivity  case  are  shown  in  Fig.  1 
for  M=50.  Note  the  thin  Hartmann  layer  near  y=H  and  the  general  inflected 
profile  in  the  core  of  the  flow.  This  inflection  is  due  to  the  rapid 
decrease  in  electrical  conductivity  with  increasing  distance  from  the 
wall  as  predicted  by  Eq.(2)  and  (3)  and  is  not  present  in  the  classical 
Hartmann  solution  where  the  conductivity  is  constant.  Such  inflected  velocity 
profiles  are  known  to  be  unstable  from  hydrodynamic  stability  considerations 
and  may  therefore  produce  a breakdown  and  subsequent  cross-flow  mixing  under 
variable  conductivity  conditions.  The  functional  form  of  the  velocity 
distribution  is  seen  to  remain  unchanged  with  varying  external  load,  a 
fact  well  known  for  constant  conductivity  calculations.  Note  also  that  the 
velocity  variation  in  the  core  goes  approximately  as  the  inverse  of  the 
electrical  conductivity  as  can  also  be  directly  surmised  from  Fq . ( 1 ) , upon 
neglecting  the  viscous  term,  which  has  significant  values  only  in  the  wall 
near  region.  That  is,  outside  the  Hartmann  boundary  layer  the  velocity 
profile  is  given  by  an  algebraic  expression  involving  a,  B„,  E and  the 
pressure  gradient  and  the  current  density  there  is  approximately  constant. 

Some  plots  of  these  analytically  determined  velocity  profiles  in  the  core 
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Details  of  the  Velocity  and  Current  Profiles  within  the  Hartmann  Layer 
for  M=1000.  The  Non-Conducting  Wall  is  at  y=H . 
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have  been  generated  for  the  different  void  distributions  given  in  Ref.[2j 
at  stations  1 through  4.  The  velocity  profiles  all  show  a characteristic 
peak  at  the  channel  center  where  the  void  fraction  is  largest. 

In  figure  2 we  show  the  velocity  and  electric  current  density  in 

our  one-dimensional  channel  flow  for  a Hartmann  number  of  M-1000  and  a 

load  factor  K=R/R  =20/21.  Again  one  notes  the  inflected  velocity  profile 

total 

and  the  extremely  thin  boundary  layer.  The  current  density  in  the  core 
of  the  flow  is  essentially  constant.  An  expanded  version  of  this  same 
solution  for  the  Hartmann  boundary  layer  is  found  in  Fig.  3 and  shows  the 
back  flow  current  in  a layer  of  thickness  Ay  = 0.0025H.  We  estimate  that 
the  typical  Hartmann  number  in  the  Argonne  liquid  metal-gas  MHD  generator [2] 
is  M=400,  so  that  the  viscous  effects  will  be  confined  almost  entirely  to 
a very  thin  layer  next  to  the  insulating  wall.  We  have  not  carried  out 
calculations  for  other  variable  electrical  conductivity  profiles  but  can 
easily  do  so  with  the  existing  computer  program.  CPU  time  for  the  velocity 
determination  for  a given  conductivity  profile  amourts  to  less  than  one 
second  and  does  not  increase  appreciably  with  increasing  M. 

An  analytical  study  has  also  been  completed  for  these  one  dimensional 
Hartmann  flows  for  the  case  of  a discrete  jump  in  electrical  conductivity 
and  viscosity  coefficient  at  a given  y between  constant  but  differing  values 
on  the  two  sides  of  the  interface.  Here  one  finds  a hyperbolic  tangent 
type  variation  in  the  velocity  W(y)  across  the  interface  with  the  transition 
thickness  proportional  to  the  reciprocal  of  the  Hartmann  number.  A current 
layer  will  generally  exist  at  such  a discontinuity  in  fluid  properties. 
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III.  TWO-DIMENSIONAL  HARTMANN  FLOWS  IN  CHANNELS  OF  RECTANGULAR  CROSS-SECTIONS 


Using  the  same  homogeneous  two-phase  approximation  as  for  the  one- 
dimensional Hartmann  flow  case,  we  next  considered  laminar  Hartmann  flows 
in  channels  of  rectangular  cross-section  with  finite  aspect  ratio.  Series 
solutions  for  such  flows  in  the  constant  conductivity  limit  have  been  obtained 
by  many  investigators  including  Shercliff  [4],  Uflyand  [5],  Chang  and  Lundgren[6] 
and  others.  However  these  calculations  are  generally  not  applicable  for  the 
wall  conductivities  encountered  in  MHD  generator  configurations  nor  are 
the  series  solutions  obtained  well  suited  for  numerical  evaluation  at  high 
Hartmann  numbers.  Accordingly,  instead  of  using  the  series  solution  method 
which  is  unlikely  to  work  for  variable  conductivity  conditions  anyway,  we 
went  directly  to  a digital  solution  approach  using  a two  dimensional  version  of 
a variable  grid  spacing  such  that  the  number  of  points  within  the  boundary 
layers  of  the  channel  walls  equal  approximately  those  within  the  core  of  the 
flow.  The  governing  momentum  and  magnetic  induction  equations  governing 
such  variable  property  Hartmann  flows  were  derived  and  in  Cartesian  coordinates 
assume  the  form 
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where  B„  is  the  constant  external  magnetic  field  applied  in  the  y direction, 

W(x,y)  the  sought  after  axial  velocity  distribution  and  B the  induced  field 

in  the  z direction.  Here  P represents  the  total  pressure  and  n=l/p  a is 

e 

the  magnetic  viscosity  related  to  the  magnetic  permeability  and  electrical 
conductivity.  We  note  that  Eqs.  (5)  and  (6)  reduce  respectively  to  the  classical 


in 


MHD  momentum  and  induction  equations [3]  in  the  limit  of  constant  h and  n. 

The  important  non-dimensional  parameters  in  the  problem  are  the  aspect 
ratio  r=b/a  of  the  channel,  where  2b  is  the  distance  between  thin  finite 
conductivity  electrodes  at  x=+b  and  2a  the  distance  between  the  two 
insulating  sides  of  the  channel  at  y=+a,  and  the  Hartmann  number  here 
defined  at  M=aB0^  o0/ho where  y0  and  o0  are  the  fluid  viscosity  and 
electrical  conductivity  at  the  wall,  respectively.  The  boundary  conditions 
appropriate  to  an  MHD  generator  channel  and  those  used  in  our  calculations 
were  that  W vanishes  on  the  channel  walls,  that  B=Const.  at  the  insulating 
walls  and  that  both  the  tangential  component  of  the  electric  field  and  the 
normal  component  of  the  magnetic  intensity  be  continuous  at  the  electrodes. 
For  the  finite  thickness  and  finite  electrical  conductivity  electrodes  under 
consideration  this  leads  to  a mixed  boundary  condition  at  the  electrodes 
identical  to  that  used  by  Chang  and  Lundgren[6],  An  important  parameter 
entering  the  boundary  conditions  at  the  electrodes  is  the  non-dimensional 
electrode  resistivity,  K^=resistivity  of  electrode/resistivity  of  fluid  plus 
electrode . 

Numerical  calculations  using  a variable  grid  spacing  based  upon  a 
boundary  layer  thickness  of  a/M  at  the  insulating  surfaces  and  one  of 
a/M  at  the  electrodes  was  carried  out  using  the  extended  Peaceman-Rachford 
alternate  direction  implicit  method [7).  The  method  involves  a relaxation 
scheme  in  which  a time  derivative  is  introduced  into  the  equations  and  the 
solutions  relaxed  until  all  time  dependence  vanishes.  The  advantage  of  this 
method  is  that  all  iterations  are  solved  implicitly,  considering  all  nodes 
simultaneously.  This  gives  numerical  stability  to  the  calculations  and  the 


necessary  matrix  inversion  is  reduced  to  inverting  a tridiagionial  system. 


Because  of  the  symmetries  inherent  in  the  problem  our  calculations  were  con- 
fined to  only  the  first  quadrant  of  the  rectangular  cross-section.  The 
velocity  and  magnetic  intensity  was  evaluated  at  A00  mesh  points  within 
this  quadrant  and  typicaly  involved  six  seconds  of  computer  time  at  a Hartmann 
number  of  20.  Typically  35  to  50  iterations  were  required  for  convergence. 

<:ur  first  calculations  were  done  for  the  constant  electrical 
conductivity^ and  viscosity  case  at  an  aspect  ratio  r=b/a=2,  a Hartmann  number 
20,  zero  load  factor  and  electrodes  of  infinite  conductivity  (K^=0) . Results 
of  the  calculation  are  shown  in  Figures  A and  5.  The  first  of  these 
shows  contours  of  constant  velocity  in  increments  of  one-eighth  the  maximum 
velocity  at  the  channel  center.  Note  the  relatively  thick  boundary  layer 
at  the  electrodes  compared  to  the  narrow  Hartmann  layer  at  the  insulating 
walls  which  are  perpendicular  to  the  applied  magnetic  field.  The  induced 
current  lines  (which  are  proportional  to  B)  are  seen  to  lie  essentially 
parallel  to  the  non-conducting  walls  with  no  appreciable  return  current  in 
the  Hartmann  layer.  This  fact  stems  from  the  zero  external  load  considered 
in  this  particular  calculation.  Large  return  currents  near  the  insulating 
walls  were  found  for  large  load  factors  as  demonstrated  in  Figs.  6 and  7. 

Note  in  Fig.  7 that  all  current  lines  remain  in  the  liquid  as  the  electrodes 
are  non-conducting  (K^=l).  To  our  knowledge  these  calculations  represent 
the  first  finite  difference  solutions  for  the  Hartmann  generator  geometry 
and  should  be  extendable  without  major  problems  to  other  aspect  ratios  and 
larger  Hartmann  numbers.  As  In  the  one  dimensional  case,  the  current  density 
within  the  core  of  the  flow  is  essentially  constant. 

Some  calculations  were  also  done  for  a variable  conductivity  fluid 
of  assumed  constant  viscosity.  For  a choice  of  conductivity  we  considered 
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Fig.  4:  Velocity  Contours  for  Two-Dimensional  Hartmann  Flow  in  the  MHD  Generato 

Configuration.  Constant  Electrical  Conductivity  and  Zero  External  Load 
Aspect  Ratio  r=b/a=2.  Hartmann  Number  M=20. 


INSULATING  WALL 


Fig.  5:  Current  Streamlines  for  Hartmann  Flow  in  the  MHD  Generator  Configuration. 

Note  the  Absence  of  Return  Currents  in  the  Hartmann  Boundary  Layer 
for  the  Zero  Load  Factor  Case  Considered. 
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Fig.  7:  Current  Streamlines  for  Hartmann  Flow  with  a Load  Factor 
and  Non-Conducting  Electrodes  (K  =1.0). 


the  hyperbolic  tangent  profiles 

a = 0o/ { 1- ( 1/4) (1-A) [1+tanh  8(l-x/a) ] [1+tanh  8(0.5-y/a)]}  , (7) 

where  A Is  an  adjustable  parameter.  Such  a conductivity  profile  is  an 
approximation  to  what  can  be  expected  in  a two-phase  MHD  channel  flow  i'n 
which  the  conductivity  near  the  wall  is  large  but  becomes  quite  small 
near  the  axis  of  the  channel  where  most  of  the  gas  will  be  located.  A 
calculation  for  A=10  is  shown  in  Figs.  8 and  9.  The  load  factor  here  is 
K=0.95  and  the  Hartmann  number  remains  at  M=20.  Comparing  these  results 
with  those  obtained  for  constant  electrical  conductivity  indicates  the 
presence  of  a high  velocity  core  near  the  channel  axis  followed  by  a velocity 
plateau  at  intermediate  distances  from  the  wall  and  then  the  usual  rapid 
velocity  decrease  in  the  boundary  layers.  Note  that  this  time  there  is  a 
large  current  backflow  in  the  vicinity  of  the  insulating  walls.  This  is  the 
result  mainly  of  the  large  load  factor  considered. 

We  have  not  at  this  point  extended  our  calculations  to  larger  Hart- 
mann numbers  in  these  two-dimensional  calculations  because  of  the  increased 
computer  time  required.  Ways  to  reduce  the  CPU  time  for  larger  M are  now 
under  investigation  including  the  possibility  of  a boundary  layer  hookup 
with  an  inviscid  core  calculation.  The  calculation  method  could  be  used  to 
make  parametric  studies  with  different  viscosity  and  conductivity  distri- 
butions should  this  seem  advisable,  and  also  used  to  determine  volume  flow 
rate  as  a function  of  channel  aspect  ratio.  It  is  likely  that  the  largest 
flow  rates  for  a fixed  area  4ab  will  occur  in  channels  of  aspect 
ratio  r=b/a  less  than  unity  as  the  boundary  layers  at  the  electrodes  are 
generally  thicker  than  the  Hartmann  layers  at  the  insulating  sides. 


3QDai3313 


i 


INSULATING  WALL 


Fig.  9:  Current  Streamlines  in  MHD  Channel  Flow  for  the  Variable  F loctrical 
Conductivity  Distribution  given  by  Eq.(7)  with  A=10. 


IV.  TWO-PHASE,  ZERO-SLIP  MODEL  FOR  COMPRESSIBLE  MHD  FLOW 
In  the  two  problems  on  Hartmann  flow  described  above  it  has  been  assumed 
that  there  was  no  variation  in  conditions  in  the  streamwise  direction.  This  is 
generally  not  true  in  a two-phase  MHD  generator  of  the  type  under  consideration 
at  Argonne  Laboratories  where  the  gas  in  the  two-phase  flow  undergoes  consider- 
able expansion  between  the  entrance  and  exit  part  of  the  MHD  channel.  Accordingly, 
we  have  also  devoted  some  time  to  examining  existing  one-dimensional  two- 
phase  slip  models  for  which  computer  codes  have  been  developed  by  the  Argonne 
group  [2]  and  attempted  to  develop  some  simplified  version  of  their  model  which 
is  consistent  with  the  experimental  observations  and  allows  for  time-variations 
in  the  variables.  In  particular,  we  examined  the  slip  model  used  by  them  and 
feel  that  there  are  questions  concerning  the  applicability  of  a churn  turbu- 
lence term  which  is  dependent  on  the  presence  of  an  average  slip  between  the 
gas  and  liquid  components.  Measured  variations  in  void  fraction  in  the  cross- 
' stream  direction  could  account  for  an  apparent  slip  predicted  by  a one- 

dimensional approximation  although  locally  there  may  be  little  velocity  difference 
between  the  components.  To  heavy  a reliance  on  such  a one-dimensional  model 
could  lead  to  incorrect  conclusions  concurring  the  upscaling  of  two-phase 
MHD  generators. 

As  a first  approximation  for  the  axial  variations  in  the  velocity,  density, 
void  fraction  and  pressure  we  have  developed  a one-dimensional  isothermal, 
zero-slip  model  in  which  the  gas  is  treated  as  ideal.  Our  objective  is  to 

» 

use  the  zero-slip  model,  whose  utility  lies  in  its  mathematical  simplicity, 
to  investigate  the  time-dependent  behaviour  of  such  flows  and  look  for 
possible  oscillatory  phenomena  that  may  be  present  and  may  lead  to  fluctuations 
in  the  power  output  in  two-phase  MHD  generators.  All  flow  properties  in  the 
cross-stream  direction  have  been  replaced  by  averages.  Within  the  framework 
of  these  approximations,  the  governing  time-dependent  equations  are  essentially 
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those  of  bubbly  two-phase  flow  as  given  by  Wallis  [9]  after  a magnetic 
force  term  is  added.  They  become 

Continuity:  (pW)  = 0 

dt  dZ 


Momentum: 

State(isothermal)  : 

Mean  Density: 
Gas  Continuity: 


r3W 

9t 

p = const,  p 


P(^7  + ^ =~  oBo(E  + WBo) 


together  with  the  subsidiary  conditions 

p = apR  + (l-cOP-j^ 

lr(p  a ^ + (p  aW)  = 0 
dt  g dz  g 


(8) 

(9) 

(10) 

(ID 

(12) 


Here  all  non-subscripted  quantities  refer  to  the  averaged  values  at 

specified  axial  positions  z.  These  equations  have  been  solved  by  a finite 

difference-predictor-connector  approach[8]  using  the  conductivity  void  fraction 

relation  given  by  eqs.  (2)  and  (3).  For  one  set  of  specified  initial  conditions 

and  boundary  conditions  based  upon  constant  pressure  at  the  channel  exit  at 

z = L and  conservation  of  momentum,  mass  and  energy  during  gas  injection  at 

the  channel  entrance  at  z = 0 and  also  closely  approximating  the  conditions 

encountered  in  the  low  temperature  NaK-N^  MHD  generator  facility  at  Argonne 

Laboratories,  one  obtaines  the  time-dependent  velocity  profiles  shown  in 

2 

Fig.  10.  For  the  large  value  of  the  interaction  parameter  N=o0B0  L/pW  used 
in  this  calculation,  the  steady-state  velocity  condition  is  seen  to  be 
approached  without  any  oscillatory  behaviour.  This  was  not  always  found  to  be 
the  case  at  lower  gas  entrance  pressures  where  transient  oscillatory  behaviour 
was  sometimes  observed.  We  are  attempting  at  the  moment  to  determine  whether 
this  oscillitory  behaviour  is  due  to  numerical  instabilities  in  the  relaxation 
procedure  used  in  the  calculations  or  represents  a real  acoustic  phenomenon 
which  may  produce  the  experimentally  observed  voltage  fluctuations  noted  under 
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certain  conditions.  In  Fig.  11  we  demonstrate  the  steady-state  axially  dependent 


pressure,  density  and  void  fraction  corresponding  to  the  steady-state  velocity 
given  in  Fig.  10.  Note,  that  as  expected  for  the  constant  channel  cross- 
section  considered,  the  velocity  and  void  fraction  generally  increase  with 
increasing  z while  the  pressure  drops. 

In  connection  with  this  zero-slip  model  we  have  also  developed  an 
approximate  analytic  model  for  channels  of  variable  cross-section.  This  model, 
in  which  we  neglect  inertia  forces  and  require  the  two-phase  flow  to  remain 
sub-sonic,  indicates  that  the  velocity  W(z)  generally  increases  with  decreasing 
pressure  and  that  approximate  constant  flow  velocities  desired  in  an  MHD  generator 
can  be  obtained  by  having  an  expanding  channel  and  tailoring  the  local  load 
factor  by  electrode  segmentation.  One  of  the  shortcomings  of  the  analytic 
model  is  that  the  acoustic  velocity  is  not  known  a priori,  being  a function  the 
gas  pressure  and  the  void  fraction. 

Some  consideration  has  also  been  given  to  the  expected  sound  velocity 
and  concomitant  potential  compressible  effects  in  the  two-phase  bubbly  flows 
existing  in  liquid  metal-gas  MHD  generators.  Employing  the  standard  approximation 
to  the  sound  velocity  in  two-phase  flows  with  void  fractions  in  the  range 
0.1  <a<0.9  as  given  in  Ref.  [9]  and  using  the  ideal  gas  law  for  assumed 
isothermal  conditions  we  find  that  the  sound  velocity  in  the  mixture  becomes 


c =| ypg/[u(l-a)p1 ] 


(13) 


where  y is  the  specific  heat  ratio  of  the  gas,  p^  the  gas  pressure,  Pj  the 
liquid  density  and  a the  void  fraction.  The  sound  velocities  predicted  by 
this  formula  yield  a minimum  at  a = 1/2  and  will  take  on  values  considerably 
lower  than  those  existing  in  the  fluid  constituents  taken  separately.  Typically, 
for  a sodium-argon  mixture  at  a channel  gas  pressure  of  76  atmospheres, 
corresponding  to  the  Argonne  Laboratories  high  temperature  facility  [ 1 ] ,the  sound 
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velocity  will  be  252  met. /sec.  This  compares  to  an  approximate  minimum 
sonic  velocity  of  24  met. /sec.  in  a water-air  mixture  at  one  atmosphere 
pressure.  In  view  of  formula  (13)  most  of  the  existing  experimental  two- 
phase  MHD  generators  operate  under  subsonic  conditions.  It  may  under 
certain  circumstances  be  of  advantage  to  operate  such  generators  in  a 
supersonic  mode  as  this  may  help  smooth  out  observed  voltage  fluctuations 
and  produce  an  acoustically  quiter  flow. 
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V.  CONCLUDING  REMARKS 


During  the  first  year  of  our  contract  work  on  two-phase  MHD 
flows  we  have  treated  the  flows  as  essentially  a homogeneous  fluid  with 
spatially  dependent  properties.  This  approach  has  yielded  informative 
new  results  which  probably  will  find  direct  application  to  some  of  the 
experimental  work  conducted  in  this  area  by  M.  Petrick  of  Argonne  Laboratories 
and  by  H.  Branover  of  Ben  Gurion  University.  We  have  treated  the  flows  as 
laminar.  This  can  be  partially  justified  by  Branovers  experimental  studies[10] 
which  indicate  that  turbulent  fluctuations  are  essentially  suppressed  when 

the  ratio  of  Hartmann  number  to  Reynolds  number  exceeds  a value  of  approxi- 

-2 

mately  10  . It  should  be  possible  to  meet  this  condition  of  turbulence 

suppression  in  experimental  two-phase  MHD  generators. 

We  plan  during  the  upcoming  contract  period  to  extend  our 
investigations  on  two-dimensional  Hartmann  flows  by  making  parametric  studies 
of  the  volume  flow  rates  through  channels  of  rectangular  cross-section  for 
various  different  conductivity  profiles.  To  handle  the  larger  Hartmann 
number  cases  without  the  use  of  an  inordinate  amount  of  computer  time,  we  will 
attempt  to  develop  a boundary  layer  solution  method  in  which  the  core  of  the 
flow  is  treated  as  inviscid.  Also  the  effects  of  varying  electrode  conductivity 
and  load  factor  will  be  further  investigated. 

The  one-dimensional,  compressible,  no-slip  model  Is  to  be  examined 
further  in  order  to  determine  the  origin  of  the  observed  oscillations  at 
lower  entrance  pressure  conditions.  We  also  intend  to  look  at  tapered  channel 
geometries  since  it  is  desirable  in  actual  two-phase  MHD  generators  to  keep 
the  flow  velocity  constant  in  its  passage  through  the  generator  duct.  An 
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attempt  will  be  made  to  include  fluid  property  variations  in  the  direction 
normal  to  the  flow  in  order  to  account  for  the  backflow  currents  existing 
in  the  Hartmann  layers  . 

We  have  not  up  to  the  present  time  considered  the  dynamics  of 
individual  bubbles  in  these  two  phase  flows,  but  plan  to  do  so  during  the 
coming  year  because  of  the  obvious  importance  of  the  suppression  of  radial 
bubble  drift  for  a properly  functioning  generator.  Analytical  work  on 
bubble  drift  in  two-phase  flows  is  nearly  non-exsistent  and  presents  a 
formidable  problem.  It  might  be  that  some  type  of  flow  separators  will 
need  to  be  installed  in  the  two-phase  flow  channel  to  suppress  such  cross- 
stream bubble  drifts. 

Finally,  some  consideration  will  be  given  to  the  time  development 
of  Hartmann  velocity  profiles  at  the  channel  entrance.  This  problem  is  of 
great  practical  importance  as  the  entrance  lengths  in  typical  MHD  generators 
may  be  comparable  with  the  actual  channel  length  so  that  a steady-state 
profile  may  actually  not  exist. 
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